#mapping.py对数据进行各种检测，
# -*- coding: utf-8 -*-
import numpy as np
from netCDF4 import Dataset
from mpl_toolkits.basemap import Basemap
import matplotlib.patches as mpatches#用于自定义图例
import matplotlib.pyplot as plt
import matplotlib as mpl
import time
np.seterr(invalid='ignore')
import cv2
plt.rcParams['font.sans-serif']=['SimHei']#显示中文标签
plt.rcParams['axes.unicode_minus']=False

def testnum():
    '''
    if i>=685:
        with open('txt.txt','a',encoding='utf-8') as f:
            f.write('{0},{1},{2},{3},{4},{5},{6},{7},{8},{9},{10},{11},{12},{13},{14},{15}\n'.
                    format(i,j,latitude[i],longitude[j],DEM[i,j],SOZ[i,j],R0_46[i,j],
                           R0_51[i,j],R0_64[i,j],R0_86[i,j],R1_6[i,j],R2_3[i,j],
                           BT3_9[i,j],BT7_3[i,j],BT11_2[i,j],BT12_3[i,j]))
            f.close()
    '''
    with open('txt.txt','a',encoding='utf-8') as f:
        f.write('{0},{1},{2},{3},{4},{5},{6},{7},{8},{9},{10},{11},{12},{13},{14},{15}\n'.
                format(i,j,latitude[i],longitude[j],DEM[i,j],SOZ[i,j],R0_46[i,j],
                        R0_51[i,j],R0_64[i,j],R0_86[i,j],R1_6[i,j],R2_3[i,j],
                        BT3_9[i,j],BT7_3[i,j],BT11_2[i,j],BT12_3[i,j]))
        f.close()
            
def getargs():
    import _test
    args=_test.degree2index(_test.setargs())
    longitude=args[0][0]
    latitude=args[0][1]
    index=args[1]

    longitude=longitude[index[0]:(index[0]+index[1])]
    latitude=latitude[index[2]:(index[2]+index[3])]
    
    filename=args[3]
    file_path=args[2][1]

    Hight,Land_SeaMask=_test.gethight(index[0],index[1],index[2],index[3])
    return longitude,latitude,Hight,file_path,filename,Land_SeaMask

#直方图均衡化
def equalizations(img):
    (b,g,r)=cv2.split(img)
    b=cv2.equalizeHist(b.astype('uint8'))
    g=cv2.equalizeHist(g.astype('uint8'))
    r=cv2.equalizeHist(r.astype('uint8'))
    return cv2.merge((b,g,r))


#调节亮度和对比度
def contrast_and_brightness(img,c,b):
    blank=np.zeros(img.shape,img.dtype)
    dst=cv2.addWeighted(img,c,blank,1-c,b)
    return dst

#G1测试
def group1(file_path,img,DEM,Land_SeaMask,flag):
    nc=Dataset(file_path)
    data=img
  
    def _R0_64():
        if SOZ[i,j]<=72:
            if (R0_64[i,j]>0.24 and dem<200) or (R0_64[i,j]>0.2 and dem>200):
                return True
            else:
                return False
        else:
            if R0_64[i,j]>(-0.014*SOZ[i,j]+1.267-0.06):
                return True
            else:
                return False

    def _R2_3():
        if dem<3000:
            if SOZ[i,j]<=56:
                if R2_3[i,j]>-0.008*SOZ[i,j]+0.676-0.06:
                    return True
                else:
                    return False
            else:
                if R2_3[i,j]>-0.008*SOZ[i,j]+0.676-0.01:
                    return True
                else:
                    return False
        else:
            return False        

    def _BT11_2():
        if BT11_2[i,j]<270:
            return True
        else:
            return False


    def BT11_2_BT_12_3():
        #   x-->sec(soz)
        #   y-->BT11.2
        def Fx22(soz,bt):
            x=1/np.cos(np.radians(soz))
            y=bt
            fxy=390.2-14.86*x-2.859*y+0.1524*x**2+0.05602*x*y+0.005239*y**2
            return fxy

        def Fx32(soz,bt):
            x=1/np.cos(np.radians(soz))
            y=bt
            fxy=238.5+107.9*x-1.895*y-21.01*x**2-0.5916*x*y+0.003557*y**2+4.524*x**3 +0.002808*x**2*y+0.001121*x*y**2
            return fxy
        
        if BT11_2[i,j]-BT12_3[i,j]>Fx32(SOZ[i,j],BT11_2[i,j]) and R0_86[i,j]<0.3:
        #if BT11_2[i,j]-BT12_3[i,j]>Fx32(SOZ[i,j],BT11_2[i,j]):
        #if BT11_2[i,j]-BT12_3[i,j]>2.7:
            return True
        else:
            return False


    def BT11_2_BT3_9():
        if dem<500 and BT11_2[i,j]-BT3_9[i,j]<-10.5:
            return True
        else:
            return False

    #沙尘检测部分
    #This is a Test for dust detection over land.
    def dust_land(dem):
        '''
        D=1.6-0.86
        d=1.38-0.86
        R1_38=(D-d)/D*R0_86[i,j]+d/D*R1_6[i,j]
        #print(R1_38)
        '''
        #数据验证部分
        def data_test():
            if all([R0_46[i,j]>0,R0_64[i,j]>0,R0_86[i,j]>0]):
                if all([BT3_9[i,j]>0,BT11_2[i,j]>0,BT12_3[i,j]>0]):
                    return True
                else:
                    return False
            else:
                return False
            
        #陆云、冰雪检测部分    
        def land_cloud():
            if any([R0_46[i,j]>0.28,_R2_3(),BT11_2_BT_12_3()]):
                return True
            else:
                return False
            
        #沙尘检测dust1
        def dust_test1():
            NDVI=(R0_86[i,j]-R0_64[i,j])/(R0_86[i,j]+R0_64[i,j])
            MNDVI=NDVI**2/R0_64[i,j]**2
            #if all([BT11_2[i,j]-BT12_3[i,j]<=-0.2,BT3_9[i,j]-BT11_2[i,j]>=15,R1_38<0.35]):
            if all([BT11_2[i,j]-BT12_3[i,j]<=-0.2,BT3_9[i,j]-BT11_2[i,j]>=15]):
                Rat1=(R0_64[i,j]-R0_46[i,j])/(R0_46[i,j]+R0_64[i,j])
                Rat2=Rat1*Rat1/R0_64[i,j]/R0_64[i,j]

                if MNDVI<0.08 and Rat2>0.005:
                    img[i,j,:]=[255,0,0]
                    return True
                else:
                    if BT3_9[i,j]-BT11_2[i,j]>20:
                        img[i,j,:]=[255,0,0]
                        return True
                    else:
                        return False
            else:
                return False
            
        #沙尘检测dust2
        def dust_test2():
            if dem<300:
                if BT12_3[i,j]-BT10_4[i,j]>=0.0010:
                    img[i,j,:]=[255,0,0]
                    return True
                else:
                    return False
            else:
                if BT12_3[i,j]-BT10_4[i,j]>=0.0100:
                    img[i,j,:]=[255,0,0]
                    return True
                else:
                    return False
 
        if Land_SeaMask[i,j]==1:
            if data_test():
                if dust_test1():
                    pass
                else:
                    if land_cloud():
                        pass
                    else:
                        if dust_test2():
                            pass
                        else:
                            group2(dem)

    #寻找地面的最大和最小高程
    def get_extreme_Height(maxh,minh):
        if Land_SeaMask[i,j]==1:
            if dem>=maxh:
                maxh=dem
            if dem<minh:
                minh=dem
        print('陆地最高海拔:  ',float(maxh))
        print('陆地最低海拔:  ',float(minh))
        return maxh,minh


    #根据海拔分级着色
    def dwHeight():
        if Land_SeaMask[i,j]==1:
            if dem<50:
                img[i,j,:]=colors[1]*255
            elif dem>=50 and dem<150:
                img[i,j,:]=colors[2]*255
            elif dem>=150 and dem<250:
                img[i,j,:]=colors[3]*255
            elif dem>=250 and dem<500:
                img[i,j,:]=colors[4]*255
            elif dem>=500 and dem<1500:
                img[i,j,:]=colors[5]*255
            else:
                img[i,j,:]=colors[6]*255
        else:
            img[i,j,:]=colors[0]*255
        return img

                 
    sizex=len(latitude)
    sizey=len(longitude)

    global i,j
    maxh=100.0
    minh=100.0
    for i in range(0,sizex):
        for j in range(0,sizey):
            check1=False
            dem=float(DEM[i,j])
            print('\r共计{}个像素,正在识别第{}个像素...'.format(sizex*sizey,i*sizey+j+1),end='')

            '''
            #子项检测
            if Land_SeaMask[i,j]==1:
                if SOZ[i,j]<=81:
                    group2(dem)
            else:
                img[i,j,:]=[255,255,255]
            '''

            '''
            if SOZ[i,j]<=81:
                #沙尘检测
                #云检测部分
                dust_land(dem)
            '''
            #HCMH测试
            if SOZ[i,j]<=81:
                if flag==0:
                    if any([_R0_64(),_R2_3(),_BT11_2(),BT11_2_BT3_9(),BT11_2[i,j]-BT12_3[i,j]>2.7]):
                        check1=True
                        
                    if check1:
                        group3(dem)
                    else:
                        group2(dem)
                else:
                    if any([_R0_64(),_R2_3(),_BT11_2(),BT11_2_BT3_9(),BT11_2_BT_12_3()]):
                        check1=True
                        
                    if check1:
                        group3(dem)
                    else:
                        group2(dem)
  
    return img

#G2测试
def group2(dem):
    check2=False

    def _R0_86_R0_64():
        if (R0_86[i,j]/R0_64[i,j])<0.8:
            return True
        else:
            if (R0_86[i,j]/R0_64[i,j])>1.6 and dem>150:
                return True
            else:
                return False

   

    def NDVI():
        NDVI=(R0_86[i,j]-R0_64[i,j])/(R0_86[i,j]+R0_64[i,j])
        if NDVI<-0.18 or NDVI>0.4:
            return True
        else:
            if dem>150:
                if (NDVI>-0.18 and NDVI<-0.14) or (NDVI>0.24 and NDVI<0.4):
                    return True
                else:
                    return False
            else:
                return False
        
    def _R0_86_R1_6():
        if dem>150 and R0_86[i,j]/R1_6[i,j]<0.75:
            return True
        else:
            return False
        
    if any([_R0_86_R0_64(),NDVI(),_R0_86_R1_6()]):
        check2=True

    '''单项测试  Cyan青色  Cyan=[0,255,255]'''

    if check2:
        #This is a clear pixel
        #do something
        #print('通过Group2,晴空区',end='')
        #img[i,j,:]=[0,255,255]
        pass

    else:
        #This is a haze pixel
        #do something
        #print('未通过Group2,雾霾',end='')
        img[i,j,:]=[255,165,0]
        #testnum()
        #pass

#G3测试
def group3(dem):
    check3=False

    def _R0_86_R1_6():
        if dem>150 and R0_86[i,j]/R1_6[i,j]<0.92:
            return True
        else:
            return False

    def NDSI():	
        NDSI=(R0_51[i,j]-R1_6[i,j])/(R1_6[i,j]+R0_51[i,j])
        if NDSI>0.36:
            return True
        else:
            return False

 
    def BT7_3_BT11_2():
        if BT7_3[i,j]-BT11_2[i,j]<-5:
            return True
        else:
            return False

    def BT11_2_BT3_9():
        if BT11_2[i,j]-BT3_9[i,j]>-9.4:
            return True
        else:
            return False

    def _BT11_2():
        if BT11_2[i,j]>288:
            return True
        else:
            return False
        
    if any([_R0_86_R1_6(),NDSI(),BT11_2_BT3_9()]):
    #if any([_R0_86_R1_6(),BT11_2_BT3_9()]):
        check3=True

    if check3:
        #This is a clear pixel
        #do something
        #print('通过Group3,晴空区',end='')
        pass
        
    else:
        #This is a cloudy pixel
        #do something
        #print('未通过Group3,云区',end='')
        img[i,j,:]=[0,255,255]

    
if __name__=="__main__":
    longitude,latitude,DEM,file_path,filename,Land_SeaMask=getargs()
    '''
    file_paths=['202103\\28\\NC_H08_20210328_0400_R21_FLDK.06001_06001_cut.nc',
                '202103\\28\\NC_H08_20210328_0400_R21_FLDK.06001_06001_cut.nc']
    '''
    file_paths=['202012\\27\\NC_H08_20201227_0100_R21_FLDK.06001_06001_cut.nc',            
           '202012\\27\\NC_H08_20201227_0300_R21_FLDK.06001_06001_cut.nc',            
           '202012\\27\\NC_H08_20201227_0500_R21_FLDK.06001_06001_cut.nc',            
           '202012\\27\\NC_H08_20201227_0700_R21_FLDK.06001_06001_cut.nc']
           
    names=['LT9','LT11','LT13','LT15']
    txts=[r'$\bf{(a)}$',r'$\bf{(e)}$',r'$\bf{(i)}$',
          r'$\bf{(b)}$',r'$\bf{(f)}$',r'$\bf{(j)}$',
          r'$\bf{(c)}$',r'$\bf{(g)}$',r'$\bf{(k)}$',
          r'$\bf{(d)}$',r'$\bf{(h)}$',r'$\bf{(l)}$']
    
    for index in range(len(file_paths)):
        #i=13
        filename=file_paths[index]
        start_time=time.time()
        name=names[index]
        #file=(file_path+filename[0]).split('.nc')[0]+'_cut.nc'
        nc=Dataset(filename)
        lon,lat=np.meshgrid(longitude,latitude)
        
        m=Basemap(projection='cyl',urcrnrlat=max(latitude),llcrnrlat=min(latitude),
                  llcrnrlon=min(longitude),urcrnrlon=max(longitude),ellps='WGS84')

        x,y=m(lon,lat)#当设置inverse=True时翻转坐标,效果相反
        extent=(x.min(),x.max(),y.min(),y.max())
        SOZ=nc.variables['SOZ']
        R0_46=nc.variables['albedo_01'][:]
        R0_51=nc.variables['albedo_02'][:]
        R0_64=nc.variables['albedo_03'][:]
        R0_86=nc.variables['albedo_04'][:]
        R1_6=nc.variables['albedo_05'][:]
        R2_3=nc.variables['albedo_06'][:]
        BT3_9=nc.variables['tbb_07'][:]
        BT7_3=nc.variables['tbb_10'][:]
        BT8_6=nc.variables['tbb_11'][:]
        BT11_2=nc.variables['tbb_14'][:]
        BT12_3=nc.variables['tbb_15'][:]
        BT10_4=nc.variables['tbb_13'][:]
        img=np.zeros((len(latitude),len(longitude),3))
        img[:,:,2]=R0_46/np.max(R0_46)*255
        img[:,:,1]=R0_51/np.max(R0_51)*255
        img[:,:,0]=R0_64/np.max(R0_64)*255
        img=equalizations(img)#直方图均衡化
        img=contrast_and_brightness(img,1.1,20)#对比度/亮度


        fig=plt.subplot(3,4,int(index+1))
        fig.imshow(img,extent=extent,alpha=1)
        fig.text(0.79,0.02,txts[3*index],transform=fig.transAxes,color='blue')
        plt.title('RGB Image '+name+':00',color='blue')
        if index==0:
            plt.ylabel('Latitude($\degree$)')

        img=group1(filename,img,DEM,Land_SeaMask,0)
        end_time=time.time()
        print('识别耗时{}s\n开始绘图'.format(end_time-start_time))

        fig=plt.subplot(3,4,int(index+5))
        plt.imshow(img,extent=extent,alpha=1)
        if index==0:
            plt.ylabel('Latitude($\degree$)')
        fig.text(0.79,0.02,txts[3*index+1],transform=fig.transAxes,color='blue')
        plt.title('HCHM算法'+name+':00',color='blue')
        cv2.imwrite('Raw'+name+'.jpg',img[:,:,::-1], [int(cv2.IMWRITE_JPEG_QUALITY),100])


        img=np.zeros((len(latitude),len(longitude),3))
        img[:,:,2]=R0_46/np.max(R0_46)*255
        img[:,:,1]=R0_51/np.max(R0_51)*255
        img[:,:,0]=R0_64/np.max(R0_64)*255
        img=equalizations(img)#直方图均衡化
        img=contrast_and_brightness(img,1.1,20)#对比度/亮度
        img=group1(filename,img,DEM,Land_SeaMask,1)
        
        fig=plt.subplot(3,4,int(index+9))
        plt.imshow(img,extent=extent,alpha=1)
        if index==0:
            plt.ylabel('Latitude($\degree$)')
        fig.text(0.79,0.02,txts[3*index+2],transform=fig.transAxes,color='blue')
        plt.title('改进算法'+name+':00',color='blue')
        cv2.imwrite('New'+name+'.jpg',img[:,:,::-1], [int(cv2.IMWRITE_JPEG_QUALITY),100])

        
    plt.figtext(0.45,0.04,'Longitude($\degree$)')
    plt.show()
